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Abstract 

We present a parameter-free Wilson-type lattice Dirac operator with an 81-point stencil 
for the covariant derivative and the Laplacian which attempts to minimize the breaking 
of rotational symmetry near the boundary of the Brillouin zone. The usefulness of this 
"Brillouin operator" in practical applications is explored by studying the scaling of pseudo- 
scalar decay constants in quenched QCD, with rather good results in the physical charm 
region. We also investigate the suitability of this operator as a kernel to the overlap 
procedure. Here, the resulting overlap operator is found to be cheaper to construct and 
significantly better localized than the variety with the standard Wilson kernel. 



1 Introduction 

Apart from being formally correct, a good lattice action should satisfy several requirements: 
(i) it should induce small cut-off effects in on-shell quantities, (ii) it should have a continuum- 
like dispersion relation, and (Hi) it should be cheap to simulate. Unfortunately, at least in 
the fermion sector these requirements tend to be in conflict with each other. For instance, the 
classic Wilson action [TJ|2] is good on (in), but not so much on the first two points. By contrast, 
an overlap action J3JH] with this operator as a kernel is significantly better on (i), but worse on 
the remaining two points. 

In the literature, there are two main avenues for obtaining a better fermion discretization. 
The "bottom-up" approach is to expand physical quantities in powers of the lattice spac- 
ing a, and to demand that the leading cut-off effects are proportional to a n a (or even a 2 ), 
where a is the strong coupling constant and n some power. This program of perturbative or 
non-perturbative 0(a)-improvement has been carried out successfully [5H9]. The "top-down" 
approach starts from the concept of a perfect action with zero cut-off effects, perfect chiral 
symmetry in the sense of the Ginsparg- Wilson relation [TO], and a continuum-like dispersion 
relation. To realize these goals exactly, a Dirac operator D(x,y) is needed with non-zero entry 
for each (x, y) pair, which is in strong conflict with criterion (Hi) above. 

In practice, one would like to maintain some degree of sparsity, that is to have an operator 
which is zero whenever x and y are further apart than a certain threshold. In the literature 
the most prominent attempts to realize ultralocal approximate derivatives of the ideal perfect 
action go by the name truncated perfect action [TTriT3] . hypercube action [T1HT6] . and chirally 
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improved action [T7JEEH]- They differ by the extent through which they make use of the full 
Dirac- Clifford algebra (the "continuum" operator uses only 7 M with fi = l..A and the identity), 
and by the criteria used to pin down the various coefficients. 

In this article we pursue a similar approach, albeit with a different focus of which properties 
should be optimized. We do not attempt to reduce 0(a 2 ) cut-off effects or the amount of 
chiral symmetry breaking, since it is known how one can get rid of these effects by adding local 
improvement terms and/or using the overlap recipe. By contrast, we strive for good overall 
appearance of the eigenvalue spectrum of D, and for a continuum-like dispersion relation, 
because for these properties no systematic improvement scheme is known. 

To ease the discussion let us consider the improved Wilson ("clover") Dirac operator 

with a^p = ^[7^, 7i/] and F^ v the hermitean clover-leaf field-strength tensor. In the wavy brackets 
there is a discrete Laplacian whose job is to lift 15 out of the 16 species, such that the resulting 
operator is doubler-free. In other words, the structure of the Wilson operator is 

D(x, y) = 53 7/^V^ td (x, y) - °- A std (x, y) + m 8 x>y + improvement term (1) 
ft 1 

where V^ 3 denotes the forward-backward symmetric covariant derivative with 2-point stencil, 
and A the standard covariant Laplacian with 9-point stencil. The mass parameters in the 
two representations above relate through 1/(2ac) = 4 + arriQ. 

The idea explored in this paper is to start from (pQ), and to replace the covariant derivative 
V^ td and the Laplacian A std by similar discretizations with improved properties. As a tribute to 
criterion (Hi) above, we shall include only one adjacent layer in each positive or negative direc- 
tion. Accordingly, both stencils have support on at most 3 d points in d space-time dimensions 
(which, in the following, will be referred to as 2D, 3D, 4D for d = 2,3,4, respectively). The 
choice of the final operator is based on a Darwinistic selection rule. Both for V M and A, a few 
varieties with distinct properties are considered, and for each combination the resulting Dirac 
operator is implemented. Based on the respective eigenvalue spectra and free field dispersion 
relations, we select the most promising combination in 2D (Sec. 2), 3D (Sec. 3) and 4D (Sec. 4). 
Fortunately, it turns out that one choice fares best regarding either criterion, and this choice is 
the same in any dimension. The resulting operator has no tunable parameters, and maintains 
the property of 75-hermiticity, i.e. 75-D75 = D*. Details of our implementation, including the 
overall link smearing strategy, the gauge covariant derivatives (based on a summation over 
all shortest paths with backprojection to the group) and tree-level clover improvement, are 
specified in Sec. 5. Practical tests in quenched QCD, with a focus on scaling studies of simple 
quantities, and in comparison to an analogously defined link-smeared tree-level clover improved 
Wilson operator, are reported in Sec. 6. In Sec. 7 we explore the suitability of our Brillouin op- 
erator as a kernel to the overlap procedure, finding a noticeable reduction of the condition 
number of the shifted hermitean kernel, and a significant improvement of the locality of the 
resulting overlap operator. A summary of our findings is given in Sec. 8, and details of all 
stencils, both in position and momentum space, are arranged in four appendices with the hope 
that they might prove useful in applications beyond lattice QCD. 
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2 Construction and main features in 2D 



2.1 Summary of 2D Laplace stencils 

The "standard" stencil of the Laplacian in 2D and the "tilted" variety (as defined in App. A) 
have the Fourier space representation (with ki = apt the dimensionless wave-number) 

a 2 A std (A;i,£; 2 ) = 2 cos(Jfei) + 2 cos(£; 2 ) - 4 

= -4sin 2 (fci/2) -4sin 2 (A; 2 /2) (2) 

a 2 A m {h,k 2 ) = 2cos(A;i)cos(A; 2 ) -2 

= 8 cos 2 (A;i/2) cos 2 (A; 2 /2) - 4 cos 2 (A; 1 /2) - 4 cos 2 (A;2/2) (3) 

respectively. From the stencil notation in App. A it is easy to see that in position space the 
"standard" Laplacian (T5]) has only 1-hop contributions (apart from the center element), while 
the "tilted" Laplacian fl3J) has only support at the edge of the 3 2 -point area around the center. 
Both of them discretize the continuum Laplacian in the sense that they deviate from the 
continuum behavior A con = —pi —p\ through 0(a 2 )-suppressed terms. Note, however, that the 
"tilted" version ([3]) differs from the "standard" variety (J2J) by having a second zero at the edge 
of the Brillouin zone, i.e. at k\ = k 2 — 71 (in the convention where the Brillouin zone ranges 
from —7i /a to n/a in every direction). 

In 2D these two stencils form a basis of all Laplace filters with (at most) a 9-point stencil. 
By taking a linear combination A(ki,k 2 ) = aA std (fc l5 k 2 ) + (1 — a) A tll (k\, k 2 ) one may try to 
improve certain properties of the discretized Laplacian. In particular, reducing the breaking of 
the rotational symmetry of the continuum operator is important. Two choices of a are popular 
in the literature. First, a = 1/2 leads to (what we call) the "Brillouin" filter [19] 

a 2 A bn (£;i, k 2 ) = cos(fci) cos(/c 2 ) + cos(/ci) + cos(k 2 ) — 3 

= 4cos 2 (A; 1 /2)cos 2 (A; 2 /2) -4 (4) 

since a 2 A hn (ki, k 2 ) = —4 whenever one of the momenta is ±71 /a. In other words, a 2 A bri takes 
a constant value on the entire boundary of the Brillouin zone. Second, the choice a = 2/3 yields 
the "isotropic" filter or stencil (see [2D] and Refs. 6-7 therein) 

a 2 A iso (k u k 2 ) = [2 cos(ifei) cos(fc 2 ) + 4 cos(An) + 4 cos(A; 2 ) - 10]/3 

= [8 cos 2 (A; 1 /2) cos 2 (A; 2 /2) + 4 cos 2 (A; 1 /2) + 4 cos 2 (A; 2 /2) - 16]/3 (5) 

since for small momenta a 2 A lso (ki,k 2 ) = — a 2 [p 2 +p 2 ] + a A [p\ +p 2 ] 2 /12 + 0(a 6 ) has 0(a 4 ) 
terms which depend only on the combination p\+p\. Put differently, the continuum relation 
^con _ _p 2 _ p 2 j g v i i a ted on-axis (to this order) in the same manner as off-axis. Note that 
this improvement strategy differs from the usual one, where one tries to remove 0(a n ) terms, 
for ever-larger n, along the axes only. 

In Fig.[TJ the momentum space representation of the four Laplacians is shown as a mesh 
plot (left) and as a contour plot (right). We choose a 24 2 lattice, and arrange the center 
of the Brillouin zone (p = 0) in the center of the frame, i.e. the boundaries correspond to 
momenta p = ±7r/a. The standard Laplacian that appears in the Wilson operator has a zero 
in the center, and decreases quadratically as one moves away from this point. As we follow 
the boundary of the Brillouin zone, it oscillates between —4 and —8. The tilted Laplacian is 
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2D: lap_std for L=24 2D: lap_std for L=24 




Figure 1: Fourier transformation of the four Laplace stencils considered in 2D. 
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rather different, since it shows a second zero at p= (n/a, vr/a), a quarter of which is seen in 
each corner. The Brillouin Laplacian has just one zero and achieves complete flatness at the 
boundary of the Brillouin zone. Finally, the isotropic Laplacian achieves best isotropy near the 
center of the Brillouin zone, meaning that one can move relatively far out from the center until 
its equipotential lines become noticeably non-circular. 

In Fig.[2]the momentum space representations of the three derivatives specified in App. A are 
shown as a mesh plot (left) and as a contour plot (right). The standard derivative that appears 
in the Wilson operator is a pure sin(pi/a), without any structure in the transverse direction. 
The Brillouin derivative modulates the transverse direction to the point that a strict zero is 
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realized on the entire boundary. The isotropic operator modulates the transverse direction in 
a less pronounced manner. For the reasons behind the name of this latter operator, which may 
sound a bit paradoxical, see App. D. 

2.2 Eigenvalue spectra in 2D 

Given the four choices of A discussed above and the three choices of V M , we can construct 12 
Dirac operators and study their eigenvalue spectra. As the gauge group is irrelevant in this 
step, we prepare a thermalized background in the U(l) gauge theory with L/a = 24 at (3 = 3.3. 

In Fig. [3] the eigenvalue spectra of the 12 operators without improvement (csw — 0) are 
shown. The Laplacian features as the row index of the panel, and the derivative as the column 
index. Out of these 12 constructions, 9 are undoubled fermion operators, while 3 yield two 
species in the continuum limit. Let us discuss the undoubled operators first. The three operators 
with A std have three branches, the left-most physical branch with the correct sensitivity to the 
topological charge of the gauge background, a doubly populated branch of wrong- chirality 
doublers near Re(z) = 2 and another species with the correct chirality near Re(z) = 4. Here the 
choice of derivative affects the spreading of the unphysical branches in the imaginary direction, 
but it leaves the topological properties of the spectrum unaffected. The three operators with 
A bn have only two branches, the left (physical) one is undoubled with the correct chirality, the 
right one includes three species, two with the wrong chirality and one with the correct chirality. 
The three operators with A lso have spectra which resemble those in the first row, except that the 
lifting of the last branch is reduced - in perfect agreement with what one expects on the basis 
of Fig.[TJ The field-theoretically most interesting spectra are those of the operators with A* 11 . 
Naively, one would expect that they yield a legal 2-flavor operator (in 2D), as the second row 
of Fig.[T] shows that this Laplacian has two zeros. With the naive derivative operator employed, 
this expectation happens to be correct; the resulting operator has two equal chirality specie^ 
which survive in the continuum limit and two doublers (again with equal, but this time wrong 
chirality) which decouple in the continuum limit. With any of the two remaining derivatives 
employed, things are a bit more involved, as the "thorn" or the "bump" in the middle or right 
panel of the second row illustrate. The point is that there is an interference^] between the 
dimension 5 Laplacian and the dimension 4 derivative; the "cross talk" phenomena in the 2nd 
and 3rd column of the second row exemplify that this may affect the structural properties of 
the fermion operator. In either case one of the would-be physical species fails to cling nicely on 
the imaginary axis for small momenta. While the operator in the 2nd column is clearly not a 
legal discretization (the "thorn" violates the property D~ip At 7 M ), the version in the 3rd column 
may represent a legal 2-flavor discretization of the Dirac operator (though, likely, with terrible 
cut-off effects). In summary, from this first overview it appears that the four operators towards 
the lower right corner of this figure seem most promising. 

In Fig. H] the same survey is repeated with tree-level improvement (csw — 1> cf. Sec. 5.3). 
Relative to the previous figure, changes seem to be mild. However, an interesting point is 
that the clover term shifts correct-type chirality branches (slightly) to the left and wrong- 
kind chirality branches (slightly) to the right. As a result, for the 9 undoubled operators the 
additive mass renormalization (the offset of the physical branch at zero imaginary part) is 

-'^Note the difference to staggered fermions in 2D, where the two species have opposite chiralities. 

2 We avoid the word "mixing", because this is a phenomenon which persists in the weak coupling limit. 
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Figure 3: Eigenvalue spectra of all operators considered in 2D with csw = 0. 
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always reduced. Also our statements in the previous paragraph regarding the chiralities of the 
unphysical branches can now be checked, because they map into a prediction of the effect of the 
clover term. At this point we can probably say that (A bn , V lso ) fares best in the sense that its 
eigenvalue spectrum is closest to that of an operator satisfying the Ginsparg- Wilson relational. 

2.3 Free field dispersion relations in 2D 

As mentioned in the introduction, the free-field dispersion relation of the fermion operator is of 
utmost importance, as this is a property for which there is no systematic improvement scheme 
(apart from taking the continuum limit). With standard 7- matrix identities it follows that the 
inverse of D = £ 7^ - \ A +m is given by D' 1 = (- £ 7^ - § A + m)/([| A -m] 2 - £ V 2 ), 
where r is the Wilson parameter. Accordingly, to work out the dispersion relation we have to 
search for zeros of [|A — m} 2 — £ V 2 , where A and V denote any one of the Laplacians or 
derivatives introduced above. 

In Fig. [5] we show, for each operator, the real solutions for r = 1 and m = over half the 
Brillouin zone on a 2D lattice with L/a = 48. The dispersion relation of the standard Wilson 
operator (A std , V std ) deviates soon from the dashed line, which corresponds to the continuum 
dispersion relation; in particular towards the boundary of the Brillouin zone the distortion is 
significant. Black boxes indicate a second real solution. If sufficiently high, this is harmless, as 
this branch decouples in the continuum limit. Note that (in certain parts of the Brillouin zone) 
some operators have only complex roots. While this proves, again, irrelevant in the continuum 
limit, it is certainly not a desirable feature. Overall, it is clear that the combination (A bn , 
V lso ) fares best in the sense that its dispersion relation is closest to the one in the continuum. 

3 Construction and main features in 3D 
3.1 Summary of 3D Laplace stencils 

The "standard" stencil of the Laplacian in 3D and the "tilted" variety (as defined in App. B) 
have the Fourier space representation 

a 2 A std {k 1 ,k 2 ,k 3 ) = 2 cos(fci) + 2 cos(A; 2 ) + 2 cos(Jfe 3 ) - 6 

= -4sin 2 (A;i/2) - 4 sin 2 (A; 2 /2) - 4 sin 2 (A;3/2) (6) 
a 2 A m (ki, k 2 , k 3 ) = 2 cos(h) cos(k 2 ) cos(/c 3 ) - 2 

= 16 cos 2 (A; 1 /2) cos 2 (A;2/2) cos 2 (A;3/2) - 8 cos 2 (A; 1 /2) cos 2 (A;2/2) - ... 
+4cos 2 (A; 1 /2) + ... -4 (7) 

respectively, with the ellipses denoting cyclic permutations. From the stencil notation in App. B 
it is easy to see that the former has only 1-hop contributions, while the latter has only 3- 
hop contributions (apart from the central element). For asymptotically small momenta they 
both reduce to the continuum relation A = pf+p 2 +pl, but the "tilted" stencil has three 
additional zeros at the boundary of the Brillouin zone [A td vanishes at (k\, k 2 , k 3 ) — (71", n, 7r)/2 = 
(±7r, ±7r, ±7r)/2 with an odd number of minus signs]. 

3 The eigenvalue spectrum of such an operator is in the unit circle centered at the point 1 on the real axis (10j . 
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Figure 5: Free-field dispersion relations of all operators considered in 2D, where |p| m ax = 7r/a- 
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In 3D the discretizations of the continuum Laplacian which are analogous to (jl]) and fl^D in 
2D are no longer simple linear combinations of ([H]) and ([7]), because one could come up with a 
Laplacian which has only 2-hop contributions (apart from the center element). They read 

a 2 A bn (fci, k 2 , ks) = [cos(fei) cos(k 2 ) cos(k 3 ) + cos(fci) cos(k 2 ) + ... + cos(fci) + ... — 7]/2 

= 4cos 2 (A;i/2)cos 2 (A; 2 /2)cos 2 (A;3/2) -4 (8) 

a 2 A iso (k l7 k 2 ,k 3 ) = [cos(fci) cos(fc 2 ) cos(fc 3 ) + 3 cos(fci) cos(fc 2 ) + ... + 5cos(£;i) + ... - 25]/6 

= [4cos 2 (A;i/2)cos 2 (A; 2 /2)cos 2 (A;3/2) +4cos 2 (/ci/2)cos 2 (A;2/2) + ... - 16]/3(9) 

respectively, and their distinctive features are as follows. The "Brillouin" Laplacian (jSJ) takes 
a constant value on the entire boundary of the Brillouin zone, since a 2 A bn (fci, k 2 , k 3 ) = —4 
whenever one of the momenta is ±7r/a. On the other hand, the Laplacian (Q is called "isotropic" 
since a 2 A iso (£;i, k 2 , k 3 ) = -a 2 [k 2 + k 2 2 + kl] + a 4 [A; 2 + A; 2 + A; 2 ] 2 /12 + 0(a 6 ) has 0(a 4 ) terms which 
depend only on the combination k\+k 2 +k 2 . In other words, A lso (A;i, k 2 , k 3 ) respects rotational 
symmetry even in the leading term through which it deviates from the continuum. 

In 3D there are 3 linearly independent Laplacians with (at most) a 27-point stencil, and 
any 3 out of the 4 elements (J6j|9j) form a basis. A systematic treatment is given in [21 J . 

3.2 Eigenvalue spectra in 3D 

Like in the preceding section, with four options for A and three for V, we can construct 12 Dirac 
operators and study their eigenvalue spectra. As the gauge group is irrelevant in this step, we 
prepare a thermalized background in the U(l) gauge theory with L/a = Yl at (3 = 2.2. A point 
worth mentioning is that in odd dimensions there is an ambiguity regarding the representation 
of the 7-matrices [22]; we opt for the 4-dimensional representation (the same one that we will 
use in 4D). 

In Fig.[6]the eigenvalues of the 12 operators without improvement (csw = 0) are shown. This 
time we refrain from showing the counterpart with improvement, as the difference is (again) 
minor. Just as in the previous section, the Laplacian features as the row index of the panel, 
and the derivative as the column index. Out of these 12 constructions, 9 are undoubled fermion 
operators, while 3 yield four species in the continuum limit. The gross features of most operators 
are rather similar to those in the 2D case, which was discussed in great detail above. Perhaps 
the most significant difference is that the operators with the standard Laplacian (first row) 
have branches at Re(z) ~ 0,2, 4, 6, with multiplicities 1,3,3,1, respectively. If the standard 
Laplacian is replaced by the Brillouin Laplacian (third row) or the isotropic Laplacian (fourth 
row), the doublers are lifted more equally; in particular the former alternative arranges them 
all near Re(z) ~2. With the tilted Laplacian and the standard derivative, the 8 species arrange 
themselves in groups of 4 at Re(z) ~ and Re(z) ~ 2, respectively. As soon as the standard 
derivative is replaced by the Brillouin or isotropic variety, some of the 4 would-be-physical 
modes cross over to the unphysical side so quickly, that the resulting operator is barely usable. 
Looking at the whole figure, one would say that the combination (A bn , V lso ) fares best in the 
sense that its eigenvalue spectrum is reasonably circular. 

3.3 Free field dispersion relations in 3D 

In Fig.UJwe show, for each operator, the real solutions for r = 1 and m = over half the Brillouin 
zone on a 3D lattice with L/a = 32. The dispersion relation of the standard Wilson operator 
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Figure 6: Eigenvalue spectra of all operators considered in 3D with csw 
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Figure 7: Free-field dispersion relations of all operators considered in 3D, where |p| max = V2n/a. 
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(A std , V std ) deviates soon from the dashed line, which corresponds to the continuum dispersion 
relation; in particular towards the boundary of the Brillouin zone the distortion is significant. 
In 3D the dispersion relation is no longer a simple curve, it depends on the orientation of the 
spatial momentum. If p is chosen on axis, the 2D dispersion relation is reproduced. The latter 
ends at y/2ir/a and features as an embedding curve to the 3D dispersion relation which now 
reaches out to y/3n/a. Again, some operators admit a second real solution (open boxes) which 
decouples in the continuum, and some operators have, for certain combinations of (pi,p 2 ), only 
complex solutions. Overall, it is clear that the combination (A bn , V lso ) fares best in the sense 
that its dispersion relation is closest to the one in the continuum. 

4 Construction and main features in 4D 
4.1 Summary of 4D Laplace stencils 

The "standard" stencil of the Laplacian in 4D and the "tilted" variety (as defined in App. C) 
have the Fourier space representation 

a 2 A std (£;i, k 2 , k 3 , h) = 2 cos(Jfei) + 2 cos^) + 2 cos(A&) + 2 cos(Jfe 4 ) - 8 

= -4 sin 2 (A;i/2) - 4 sin 2 (A; 2 /2) - 4 sin 2 (fc 3 /2) - 4 sin 2 (fc 4 /2) (10) 
a 2 A ta (fci, k 2 , k 3 , fc 4 ) = 2 cos(fci) cos(fc2) cos(k 3 ) cos (£4) — 2 

= 32 cos 2 (A: 1 /2) cos 2 (A;2/2) cos 2 ^/ 2 ) cos 2 ^/ 2 ) - 16 cos 2 (A; 1 /2) cos 2 (A; 2 /2) 
cos 2 (A;3/2) - ... + 8cos 2 (A;i/2)cos 2 (A;2/2) + ... - 4 cos 2 (A;i) - ... (11) 

respectively, with the ellipses denoting cyclic permutations. From the stencil notation in App. C 
it is easy to see that the former has only 1-hop contributions, while the latter has only 4-hop 
contributions (apart from the central element). For asymptotically small momenta they both 
reproduce the continuum relation A = pl+p^P^-pl, but the "tilted" stencil has seven additional 
zeros at the boundary of the Brillouin zone [A* 11 vanishes at (ki, k 2 , k 3 , k^) — (tt, tt, tt, 7r)/2 = 
(±7r, ±7r, ±7r, ±7r)/2 with an even number of minus signs]. 

In 4D the discretizations of the continuum Laplacian which are analogous to (HI 13) read 

a 2 A bn (fci, k 2 , k 3 , k^) = [cos{k\) cos(k 2 ) cos^) cos(/c4) + cos(fci) cos(/c2) co^k^) + ... 

+ cos(/ci) cos(k 2 ) + ... + cos(A;i) + ... - 15]/4 
= 4cos 2 (A; 1 /2)cos 2 (A;2/2)cos 2 (A;3/2)cos 2 (fc4/2) -4 (12) 
a 2 A lso (ki, k 2 , k 3 , /c 4 ) = [2cos(fci) cos(A;2) cos(k 3 ) cos(/c4) ± 7cos(fci) cos(/c 2 ) cos(/c 3 ) + ... 

+20cos(A;i)cos(A;2) + ... + 25 cos(fci) + ... - 250]/54 
= [16 cos 2 (A;i/2) cos 2 (A;2/2) cos 2 (A>3/2) cos 2 (/c 4 /2) 
+20 cos 2 (A; 1 /2) cos 2 (A;2/2) cos 2 (A;3/2) + ... 

+ 16 cos 2 (A;i/2) cos 2 (A;2/2) + ... - 16 cos 2 (A;i/2) - ... - 128]/27 (13) 

respectively, and their distinctive features are as follows. The "Brillouin" Laplacian (|12p takes 
a constant value on the entire boundary of the Brillouin zone, since a 2 A bn (/ci, k 2 , k 3 , k^) = 
—4 whenever one of the momenta is ±7r/a. On the other hand, the Laplacian ( jTBl) is called 
"isotropic" since a 2 A iso (k u k 2 ,k 3 ,h) = -a 2 [kl + kl + kl + kj] + a 4 [A; 2 + A; 2 + A; 2 + A; 2 ] 2 /12 + 0(a 6 ) 
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has 0(a A ) terms which depend only on the combination kf + k^ + k^ + kf. In other words, 
A lso (A;i, &2, &3, ki) respects rotational symmetry even in the leading term through which it 
deviates from the continuum. 

In 4D there are 4 linearly independent Laplacians with (at most) an 81-point stencil, and the 
4 elements (1 1 QUI 3 p form a basis. We are unaware of any systematic treatment in the literature. 

4.2 Eigenvalue spectra in 4D 

Like in the previous two sections, with four options for A and three for V, we can construct 
12 Dirac operators and study their eigenvalue spectra. As the gauge group is irrelevant in this 
step, we prepare a thermalized background in the U(l) gauge theory with L/a = 6 at (3 = 1.1. 

In Fig.[S] the eigenvalues of the 12 operators without improvement (csw — 0) are shown. 
Again, we refrain from showing the counterpart with improvement, as the difference is marginal. 
Following the tradition of the previous sections, the Laplacian features as the row index of the 
panel, and the derivative as the column index. Out of these 12 constructions, 9 are undoubted 
fermion operators, while 3 yield eight species in the continuum limit. Again, the gross features 
of these operators are rather similar to their 2D and 3D counterparts. This time, the operators 
with the standard Laplacian (first row) have branches at Ke(z) ~0, 2, 4, 6, 8, with multiplicities 
1,4,6,4, 1, respectively, and with alternating chiralities. Replacing the standard Laplacian by 
the Brillouin Laplacian (third row) or the isotropic Laplacian (fourth row), the lifting of the 
doublers is reduced. With the tilted Laplacian and the standard derivative, the 16 species 
arrange themselves in groups of 8 at Ke(z) ~0 (correct-type chirality) and Ke(z) ~2 (wrong- 
kind chirality), respectively!. Once the standard derivative is replaced by the Brillouin or 
isotropic variety, some "cross talk" between the marginal and the irrelevant piece becomes 
apparent. Looking at the whole figure, one would say that the combination (A bri , V lso ) fares 
best in the sense that its eigenvalue spectrum is closest to that of a Ginsparg- Wilson action. 

4.3 Free field dispersion relations in 4D 

In Fig. [9] we show, for each operator, the real solutions for r = 1 and m = over half the Brillouin 
zone on a 4D lattice with L/a = 24. The dispersion relation of the standard Wilson operator 
(A std , V std ) deviates soon from the dashed line, which corresponds to the continuum dispersion 
relation, and shows large effects of anisotropy. If p is chosen on axis, the 2D dispersion relation 
is reproduced. If p is chosen as a multiple of (1, 1, 0), the 3D dispersion relation is reproduced. 
If p is chosen as a multiple of (1, 1, 1) entries at the upper border are generated; they go out to 
\/3n/a. Again, some operators admit a second real solution (open boxes) which decouples in 
the continuum, and some operators have, for certain combinations of {pi,P2,Ps), only complex 
solutions. Overall, it is clear that the combination (A bri , V lso ) fares best in the sense that its 
dispersion relation is closest to the one in the continuum. 

4 Note the difference to naive or staggered massless fermions in 4D, where both chiralities sit on top of each 
other. For the effect of non-standard staggered mass terms and the resulting eigenvalue spectra see [231124] . 
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Figure 8: Eigenvalue spectra of all operators considered in 4D with csw = 0. 
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5 Specification of operator details in 4D 



5.1 Overall smearing strategy 

Given the results in the previous three sections, the combination of "isotropic derivative" and 
"Brillouin Laplacian" seems most attractive. In other words, our preferred operator is 

D(x, y)=J2 7^ Vf (x, y) - ~ A bri (x, y) + m 6 x , y - ^ £ a^ v 8 x , y (14) 

and below we shall refer to it as the "Brillouin operator" . 

An ingredient which has proven particularly useful in the design of fermion actions with 
small cut-off effects is link-smearing, also known under the label of "fat links" [25-30J. In the 
quenched QCD tests reported below a single step of APE smearing [31] 

V^x) = Uf E (x) = P su(3) {(l-a)I+^ ]T U„(x)U^(x+i>)U}(x+fi)Ul(x)}u^x) (15) 

with a = 0.72 is applied (for the details of the backprojection to SU(3) see subsection 5.2). The 
result is used as an input in the covariant derivative and the covariant Laplacian. The latter 
operators are made gauge covariant in the simplest possible way, by summing over all shortest 
paths, with subsequent backprojection to SU(3). For instance the hyper- diagonal connections 
(4 hops) receive contributions from 24 paths, while the cubic-diagonal connections (3 hops) 
receive 6 contributions, and the square-diagonals (2 hops) just 2. 

We use the same kind of smeared gauge links V^(x) in the construction of the derivative, 
the Laplacian, and the field strength tensor. Since this change is ultralocal and modifies only 
operators of mass dimension 5 and higher, the universality class of the action is unaffected. 
Other smearing strategies are possible, e.g. only relevant pieces or only irrelevant pieces of the 
action may be smeared. However, as we are unaware of any advantage of such more complicated 
schemes, we prefer to stay with the overall smearing strategy where all "thin" links U^x) in 
( 1141) . even if within F^, are replaced by the same kind of "fat" links V^(x). 

The goal of our quenched scaling study in Sec. 6 is to confront f[T41) with the standard Wilson 
action. To compare like with like, we will use the same smearing strategy and the same kind 
of clover improvement with csw = 1 (see subsection 5.3 below) in either case. 

For completeness let us mention that, in order to simulate full QCD with a fat-link Brillouin 
or Wilson action and an HMC algorithm, one would equip either action with a smearing that 
is tailored to this purpose (e.g. "stout/EXP" [32], "n-APE" [33], "LOG" [M], "over-improved 

stout" jug). 

5.2 Details of the projection to SU(N) 

The projection of an arbitrary NxN matrix A to SU(N) is usually defined through a projection 
to U(N), followed by a projection to unit determinant. The first projection is realized as 

P U{N) {A} = A(rfA)- 1 ' 2 = (AA^y 1/2 A (16) 

where the equivalence of the two representations follows from the singular value decomposition 
A = USV ] , with unitary U,V and S > 0, resulting in P U{N) (A) = UV ] . Since A^A and AA* 
are both hermitean, either version of (Tl~6|) requires only one eigensystem. 
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The projection to unit determinant is somewhat more involved, even if we restrict the 
discussion to unitary arguments, as suggested by the above 2-step procedure. The most naive 
recipe is to divide a given U G U(N) by the N-th root of its determinant. Unfortunately, this 
is not a valid procedure, since there is a finite (non-zero) likelihood^ that the argument has 
det(£7) = — 1, which lies on the branch cut. It is thus necessary, in general, to distribute the 
phase rotation (to go from det = e lr ^ to det = l) unevenly among the iV eigenvalues. 

In our opinion a particularly compelling option for fixing this ambiguity is to notice that 
the U(N) projection defined in (jTBJ) can be understood as the result of the recipe 

P U{N) {A} = x mm ) tr{(A-X)t(A-X)} (17) 

and to define the complete projection by using this recipe for the SU(N) group, i.e. via 

Psu(N){A} = y mm N) tx{(A-Yy(A-Y)} (18) 

in a single step, where an algorithmic solution has been proposed in [36J: 

1. Perform a singular value decomposition A = U SV^ with U, V G U (N) and S > a diagonal 
matrix with positive entries. Indeed, X — UV^ is the projection to U(N), but det(X)^l. 

2. Compute det( A) =pexp(i0). Incidentally, det(S) =p and det(UV^) = exp(i0). The matrix 
U exp(—i4>/N)V* is in SU(N), but it is, in general, not the one which is closest to A. 

3. Find the solution for the phases of the matrix D = diag(exp(i6 l i), exp(i#jv))j subject 
to the constraint J2^i J r4 > = ^ (mod 2ir), which maximizes HeTi{A^U DV^). By means of 
the original singular value decomposition, the latter expression equals ReTr(SD), and the 
expression to be maximized^ is s x cos(6*i) + ... + sn cos(9 n ), still subject to the constraint 
X]#i + = O (mod 2n). The matrix Y = UDV^ is the desired solution. 

From (|18p it is clear that, if A is subject to a random gauge transformation A — > g\Ag\ with 
gip G SU(N), the effect must be that the solution Y transforms as Y — > g{Yg\- Up to a set 
of gauge configurations of measure zero, the singular value decomposition A = USV^ is unique 
(here we assume a specific ordering in S, e.g. Si > ... > s N > 0). As a result, the effect of the 
random gauge transformation is just U — > giU, V — » g 2 V, while S and the expression to be 
maximized are unchanged, and the net effect is thus UDV^ — > g\UDV^g\, as expected. 



5.3 Tree- level improvement 

The Symanzik effective field theory of cut-off effects of undoubted lattice Dirac operators is 
based on an analysis of all local mass dimension 5 operators consistent with the symmetries 

5 The issue arises if A is the sum of only two SU(N) matrices. The determinant of A = U + V is real, since 

det(U+V) = det(J7(F t + t/ t )V r ) = det^ + U^ = [det(U+V)}* 

with no generalization if A is a sum of three and more unitary matrices. Hence, with arbitrary U, V € SU(N) 
it may happen that det(U+V) lies on the negative real axis. Taking a look at (fT6|) one realizes that det(A) <0 
implies det(£>) = — 1, where B = Pu^{A} and A the original sum of two SU(N) matrices. 

6 Note that the global maximum is required. A possible strategy is to perform a scan, in regularly spaced 
intervals, over all unconstrained 9i, followed by a local maximization starting from the largest value obtained 
in the survey. Thus, staying content with intervals A#i = 7r/3, the global scan requires function calls. 
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of the theory [5HH]. As long as one is content with perturbative or non-perturbative 0(a)- 
improvement of on-shell Green's functions, contact terms can be ignored and it suffices to 
add a single improvement term, the so-called clover term which is included in (IT| H%|) . Going 
through the arguments of [5H9], one realizes that the leading contribution, in the weak coupling 
expansion, to the coefficient in front is independent of the details of the covariant derivative 
and Laplacian occurring in the operator. In other words, csw = 1 holds true, at tree level, also 
for our Brillouin operator flHJ), while subleading contributions are, of course, different. 

With csw = 1 the leading cut-off effects are 0(aa), and the overall smearing strategy does 
not change this. However, due to the smearing the coefficient of the 0(aa) term might be so 
small, that the formally subleading 0(a 2 ) cut-off effects might prove numerically dominant. 



6 Practical tests in quenched QCD 
6.1 Scale setting and overall tuning strategy 

We use the Wilson gauge action and a parametrization of r^/a consistent with asymptotics 

, , . , 4?r 2 „ 1 - 8.2384//3 + 15.310//3 2 , , 

logfoAO = 33 ! _ 2.7395//? - 11,526/^ (19) 

which is based on data from [38]. Upon choosing L/a = 10, 12, 16,20,24 and requesting that 
L/ro = 3.2653 (tantamount to L~1.6fm if ro~0.49fm is assumed), we find that we should use 
the (3 values listed in Tab. CD Note that the quenched scale ambiguity of ~5% does not limit 
our ability to match boxes in terms of L/tq. 

We aim for comparisons at a fixed value of the light and the strange quark mass. This 
will be achieved by tuning the pion and the ss mass (without disconnected contributions) to 
(r M 7r ) 2 ~1.56~1.25 2 and (r M 5s ) 2 ~ 4.56, respectively. Given that M$ 8 ~2M%-M%, this will 
correspond to (r M^) 2 ~ 3.06 ~1.75 2 , and hence to a pion mass of about 500 MeV and a kaon 
mass of about 700 MeV. Finite volume effects are expected to be small, since M 7r L~4.08. 



6.2 Determination of ft cr } t and choice of K]i g u, ^strange, ^charm 

As a first step we determine K crit for either operator (with 1 APE step and csw — 1) over a range 
of (3- values, with results given in Tab. [T] and Fig.[TDl The perturbative formula reads (f3 = 6/g 2 ) 

am crit = S = Sl°F S + °(So) [< 0] (20) 

with Cf = 4/3 and om cr i t = l/(2/c cr it) — 4. At 1-loop order one finds S = 31.98644 for thin- 
link Wilson operator with csw — 1, and S = 4.07175 for the 1APE (a = 0.72) variety with 
csw = l [30], while for the Brillouin operator no perturbative information is available. 

To see how far from the perturbative regime we are, we fit our data to the rational ansatz 

- am crit = Cl f I ± C29 J (21) 

with 2 degrees of freedom, and compare the fitted C\ to the 1-loop prediction S'/(127r 2 ), where 
available. For the Wilson operator c\ = 0.0623(82) deviates significantly from the prediction 
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(3 geom. a 1 [GeV] 


K crft' /btd ( "Standard" ) 


Ct bri ("Brillouin") 


5.72 10 3 x20 1.236 
5.80 12 3 x24 1.479 
5.95 16 3 x32 1.978 
6.08 20 3 x40 2.463 
6.20 24 3 x48 2.964 


0.134516(65) 0.134533 
0.132673(47) 0.132650 
0.130760(59) 0.130769 
0.129818(45) 0.129864 
0.129362(57) 0.129303 


0.129780(64) 0.129798 
0.128594(30) 0.128582 
0.127469(48) 0.127471 
0.126940(30) 0.126973 
0.126725(42) 0.126676 



Table 1: Summary of n CTit for either operator with 1 APE step and Csw = 1, as determined from 
direct measurements at these couplings (with error bars) and through the fit JUP- 
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Figure 10: Additive mass renormalization versus g 2 for the Wilson operator and the Brillouin 
operator (both with csw = l an d Q!APE = 0.72j. In either case a rational fit with the ansatz < f21]) 
is included. Error bars are significantly smaller than the size of the symbols. 



0.0344, while for the Brillouin operator c\ = 0.0143(54), without a perturbative prediction to 
compare to. Given the quality of these fits, the interpolated K crit are more accurate than the 
direct measurements, and this is why we include these values in Tab.fD 

In order to perform a quenched scaling study we define, for each 0, three reference K-values 
which realize (roM n ) 2 = 1.56, (vqMss) 2 = 4.56 and (r M 5c ) 2 = 46.5. We call them /tn g ht, K strange , 
and Kcharm, respectively (even though the first two are heavier than the respective physical 
flavors, and the last one is lighter than the physical charm quark). These values are determined, 
for each coupling, by interpolating the results of a few tuning runs. The three reference K-values 
are then evaluated on the full ensembles, and the resulting (r Mp) 2 are compared to the target 
values in Fig. [Til It seems the tuning is accurate enough, so that we can proceed with a study 
of the scaling of the decay constants. 
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Figure 11: Summary of the final (roM^) 2 , (r M gs ) 2 , (r M 5 c) 2 to test how accurately the target 
values 1.56, 4.56, 46.5 were reached. Error bars are smaller than the size of the symbols. 



6.3 Comparing the scaling of decay constants at fixed r^Mp 

The decay constants F n , F Ss , F 5c are determined from the improved and renormalized current 

A™ = Z A (l + b A am w )(A fl + ac A d,P) (22) 

where and P denote the naive axial-vector current and pseudo-scalar density, respectively, 
and m w = mo — m cr it . In practice, m w in fl22l) is often replaced by the PCAC quark mass 



m 



pcac = EAd4A^x) + ac A d 4 P(x)]O(0)) 
2£ x (P(x)O(0)> 



(23) 



where d denotes the symmetric derivative, and usually = P is chosen to get maximal signal. 
Here it is assumed that the two quark masses are equal; in general the improvement factor in 



( 1221) is (l + b A a(rrij +mj: )/2) for flavors j, k, and the l.h.s. of 

C4 : 



is m 



PCAC 



+m p k CAC 



)/2- 



(J23J) 

We use the tree-level improvement coefficients b A = 1, c A = 0. The 1-loop renormalization 
constant Z A = l—g 2 z A / (127r 2 ), which is needed for consistency, is known for the Wilson operator 
(z A = 2. 42423 with 1 step of ckape = 0.72 smearing and csw = l is found in [30]), but not for the 
Brillouin operator. In Fig. [12] we plot the decay constants F n , F Ss , F 5c versus aa (left) and a 2 
(right). Here everything is made dimensionless through tq. In case of the Wilson operator the 
lattice-to-continuum matching factor Z A is included, but it brings a rather small shift, since it 
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Figure 12: Decay constants F n (top), F Ss (middle), F 5c (bottom) in r units versus aa (left) 
and a 2 (right). Open symbols indicate the bare values, filled symbols include the 1-loop 



is already close to 1 in the range of couplings where we have data, and it approaches 1 as a — > 0. 
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Figure 13: The ratios F Ss / ' (top) and F 5c / F Ss (bottom) versus aa (left) and a 2 (right). The 
linear fits include only 4 lattice spacings, and favor the pure a 2 over the pure aa extrapolation. 



With hindsight we can thus anticipate that also for the Brillouin operator the data without Za 
are indicative of the approach to the continuum. Comparing the two operators without Za, we 
see little difference in the light and strange pseudo-scalar data (top and middle), while there is 
a pronounced difference in the charm sector (bottom row). Hence, for the scaling of roF 5c the 
Brillouin operator seems to bring a significant improvement. 

To get rid of the Za factors, we also consider the scaling of the ratios F Ss /F n and F 5c /Fs S , 
as shown in Fig.[13j Again, we plot the data against aa (left) and a 2 (right). For the strange- 
to-light ratio all data happen to be essentially flat, so there is no advantage of one operator 
over the other. For the charm-to-strange ratio, the situation is different. Fitting the data on 
the four finer lattices with a pure aa ansatz yields two continuum extrapolated results which 
are not consistent (lower left panel). Fitting the same data with a pure a 2 ansatz leads to two 
continuum extrapolated results which are almost consistent (lower right panel). If we restrict 
the fits to the three finest lattice spacings, the values obtained with the pure aa hypothesis 
stay inconsistent, while the continuum results with the pure a 2 hypothesis become consistent. 
To prevent any misunderstanding, let us emphasize that we think that both operators have a 
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Figure 14: Fit of the mixed aa plus a 2 ansatz ( T23j) to the ratio F Ec /F Ss with 4 (left) or 5 (right) 
lattice spacings included. 



contribution in era and a 2 at accessible lattice spacings. Still, to the best of our knowledge, 
this is the first figure which indicates that, for a tree- level improved operator with some link- 
smearing, the pure a 2 hypothesis might be closer to the truth than the (formally correct) pure 
era hypothesis. Of course, with infinitely precise data one could separate the two contributions. 
To see how far we are from this ideal world, we try a fit of the ratio F 5c /F Ss with the ansatz 

F- cc /F Ss = d + d l a(a)a + d 2 a 2 (24) 

giving results shown in Fig.[T3J The fitted d\,d% of the Brillouin operator are significantly 
smaller than those of the Wilson operator. Also by looking at the fits one would say that the 
Brillouin data alone leave little doubt that the correct continuum value is somewhere near 1.85, 
while with the Wilson data alone this is far from obvious. 



6.4 Comparing the l/neiCGstab distributions at fixed r^M^ 

In quenched QCD with Wilson fermions so-called exceptional configurations (on which the 
massive Dirac operator D m could not be inverted) hindered the approach to light quark masses. 
In full QCD the functional measure suppresses configurations on which D m has near- zero modes. 
Still, the issue persists in the form of instabilities in the HMC evolution. 

In [39] it was shown that the stability of these simulations is linked to the distribution of 
the lowest eigenvalue of Dl n D m . The latter is roughly Gaussian distributed, and the simulation 
is deemed safe as long as the center of the distribution is at least four standard deviations 
away from zero. The BMW collaboration noticed that the smallest eigenvalue of D\ n D m is 
directly related to the number of iterations in the inversion, and used the inverse iteration 
count 1/uqq in the monitoring [ID]. In Fig.[TS]we present l/n B iCGstab f° r either operator at the 
values (roM^-) 2 = 1.56 and 0.56 (M m ~ 500 MeV and 300 MeV). In either case an inversion with 
the Brillouin operator requires about 60% of the forward application^ of the Wilson operator. 

7 For fixed M„ the smallest eigenvalues of the two A — D^Dm are approximately equal, while the largest 
eigenvalue is near 2.5 2 for the Brillouin operator and near 7.5 2 for the Wilson operator. Since jiqq cx y/CN(A) 
one would expect the relative iteration count to be around 1/3 for CG and around l/v3 ~ 0.6 for BiCGstab. 
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Figure 15: Distribution of the inverse iteration count 1/riBiCGstab to reach a norm e — 10 7 of 
theresidual at/3 = 6.20 with k tuned to have (r M^) 2 = 1.56 ("500MeV") or0.56 ("300MeV"). 
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Figure 16: Distribution of the logarithm of the correlator P(T/4)P(0) at /3 = 6.20 with k tuned 
to have (roM n ) 2 = 1.56 in either case. 



Finally, according to the safety criterion mentioned above, there seems to be a slight advantage 
for the Brillouin operator at low M n . Using another /3- value did not bring any major change. 

6.5 Comparing the statistical fluctuations at fixed roM n 

Reaching the same statistics for Wilson and Brillouin data may or may not be a good guide 
to obtain equally precise physics results. To compare the fluctuations with either operator we 
compare the variance in the correlator C P p(t) at t = T/4. The result is shown in Fig. [16] in 
the form of a histogram, with either k tuned to realize (roM^) 2 = 1.56. Essentially, there is no 
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noticeable difference between the two operators. Looking at other /3-values we arrived at the 
same conclusion 

7 Suitability as overlap kernel 

7.1 Details of the overlap action 

Given any undoubted (flavor symmetry respecting) "kernel" Dirac operator D^ m at a quark 
mass m, the massless overlap operator is defined through [3] 

D ov = D ovfi = P - [1 + D^ p/a (Dl n _ p/a D^ p/a )- l l 2 } = P - [1 + 75 sign( 75J D kni _ p/a )] (25) 

with 0<p<2. Traditionally, the kernel parameter p was tuned to a value above 1 to maximize 
the locality of D ov on coarse lattices [JT] • However, on fine lattices (and with some link smearing 
or filtering of the kernel also on relatively coarse ones) maximal locality is obtained for p < 1 
[i2lPf3] . As p is part of the action definition, it is desirable to keep it fixed, and we stay with 
the canonical choice p — 1 . 

The massive overlap operator follows by adding a "chirally rotated" scalar term [H] 

(X CLTYl 

D ov , m = D ov + m(l- — A, v ) = (1- ~^) D o, + m (26) 

which yields an operator with a circular eigenvalue spectrum of radius p — am/2 around the 
point (p+am/2,0) in the complex plane. 

There is still a choice to be made regarding the filtering of the underlying kernel operator 
(we use 1 and 5 APE steps) and whether one wants to equip it with a clover term. 

7.2 Comparing the near-normality of the kernels 

Wilson-type operators are usually non-normal, i.e. [D, D'] y^O jl5]. This means that the spectral 
representation takes the form D = J2^n\ ( fin)('4 ) n\-> with no simple connection between \<p n ) and 
(ip n \ and, as a result of this, no simple connection between the eigenvalue spectra of D and 
D^D. Chiral operators are usually normal, i.e. [D, D^} = for a staggered or overlap Dirac 
operator. This means that D = J2 A n |V'n)(V , n|, with \ip n ) and (ip n \ being the complex conjugate 
transpose of each other, and the spectrum of D^D can be inferred from the one of D. 

In this sense one may understand the non-normality of a Wilson-type fermion, defined as 
the norm of the commutator, as a measure of "how far" it still is from a formulation with 
continuum-like features. Therefore, we measure ||(£)tD — DD^)i]\ \ for a few dozen normalized 
Gaussian random vectors 77, with D being the Wilson or the Brillouin operator. By doing this 
on 15 configurations for each ensemble, we obtain the data shown in Fig.[T71 Unsurprisingly, 
with either operator the non-normality decreases towards the continuum, but the Brillouin 
operator fares significantly better. Surprisingly, switching on the clover term increases the non- 
normality, but it remains true that with the Brillouin operator the norm of the commutator is 
about an order of magnitude smaller than with the Wilson operator. 
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Figure 17: Non-normality \ \ [D, D'] || versus g\, after 1 APE step, for the Wilson (left) and 
Brillouin (right) operator, with csw = (top) and 1 (bottom). Note the difference in scale. 



7.3 Comparing the Ginsparg- Wilson violation of the kernels 

The Ginsparg- Wilson relation for a massless D reads [10] 



l5 D + D75 = -D^D 
P 



(27) 



and we intend to plug in our operators with k set to K C rit- Since the latter are known only for 
csw = l (cf. Tab. [I]) we do this with improvement. We measure IK-D75 + 75-D — Dj 5 D)rj\ \ for a 
few dozen normalized Gaussian random vectors rj on 15 configurations of each ensemble. The 
result is shown in Fig.fTSJ A priori it is not clear whether it makes sense to plug a distinctly 
non-chiral operator into the Ginsparg- Wilson relation (1271) . but the result of our experiment 
seems to suggest that at least for the Brillouin operator it does. 
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Figure 18: Violation of the Ginsparg-Wilson relation ( T2Tj) with p = l as a function of 6/(3 for 
either operator (see text for details). Note the difference in scale. 



7.4 Comparing the condition number of the hermitean kernels 

The cost of the overlap construction is determined by the smallest mode (in absolute magnitude) 
of the shifted hermitean kernel H kn _ p / a = 75-D kn _ p / a = 7 5 (Z) kn — p/a). Equivalently, one can 
look at the condition number of the squared operator A = _ p / a D kn In practice, one 
considers the so-called Ritz eigenvalues, i.e. the eigenvalues of the symmetric tridiagonal matrix 
that emerges from the Lanczos process on A. They approximate the extremal eigenvalues of A. 

In the top panel of Fig. [19] we plot the smallest and the largest Ritz eigenvalue of A made 
from the standard Wilson kernel (left) or the Brillouin kernel (right) as a function of the mass 
airiQ = —p. In the bottom panel the 10th and (again) the largest eigenvalues are shown. In either 
panel, 16 configurations of our finest ensemble (/3 = 6.20) are used, after 1 step of APE smearing 
is applied, and the clover coefficient is set to zero. Note that the gap between the largest and 
the smallest eigenvalue is just the condition number of A, and the gap between the largest and 
the 10th eigenvalue is the condition number of A restricted to the subspace orthogonal to the 
lowest 9 eigenmodes. Hence, after O(10) eigenmodes are projected, the Brillouin kernel allows 
for a reduced order of the polynomial or rational representation of the sign function, since its 
spectral range is 1-2 orders of magnitude smaller. 

In Fig.[20]the same exercise is repeated with 5 steps of APE smearing. The overall picture is 
unchanged; again the resulting condition number of the Brillouin kernel is significantly smaller. 
We also find little impact of 19 (instead of 9) eigenmodes projected, and whether the kernel 
is Symanzik improved or not. In short, the reduction of the condition number comes predom- 
inantly from the lowering of the largest eigenvalue (in line with what one would expect from 
the eigenvalue spectra shown in Fig.[8]). By chance, one of the configurations used in Fig. [20] 
happens to be close (in configuration space) to a barrier between two topological sectors. With 
the Wilson kernel the crossing occurs near p=1.7, with the Brillouin kernel close to p = l.l. 

In summary, we find that the shifted Brillouin kernel has a significantly reduced condition 
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Figure 19: The 1st (top), 10th (bottom) and largest Ritz eigenvalue of D\ n D m for the Wilson 
and the Brillouin operator, on 16 configs at /3 = 6.20, with csw = and 1 step of APE smearing. 
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Figure 20: The 1st (top), 10th (bottom) and largest Ritz eigenvalue of D\ n D m for the Wilson 
and the Brillouin operator, on 16 configs at (3 = 6.20, with csw — and 5 steps of APE smearing. 
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Figure 21: Localization of the p = l overlap operator with the standard Wilson kernel (left) or 
the new Brillouin kernel (right) on a free 48 4 lattice, for four directions of the separation. 



number, in particular with a bit of link smearing and after O(10) eigenmodes are projected. 
This allows for a lower degree polynomial or rational representation of the sign function. 

7.5 Comparing the locality of the resulting overlap actions 

The locality of the overlap action with standard Wilson kernel was first studied in |41j. In [46J 
it was shown that a nearly chiral (but still ultralocal) kernel can significantly improve the 
coordinate-space locality of the resulting overlap action. In [4^||4^] it was shown that even 
a slight modification through some link-smearing can lead to a considerable improvement. 
Therefore, one may hope that trading the Wilson kernel for the Brillouin kernel leads to a 
noticeable improvement of the locality of the overlap operator. 

The localization of the overlap made from the Wilson or the Brillouin kernel is shown for 
a 48 4 lattice in the free field case in Fig.[2TJ The Frobenius norm of D(x,y) is plotted as a 
function of the Euclidean distance d,2 = \ \x—y\\2- Evidently, the Brillouin kernel diminishes the 
anisotropy effects and makes the operator fall off at about twice the rate as before. 



8 Summary 

We have introduced an ultralocal single-flavor lattice Dirac operator, based on the gauge covari- 
ant versions of V lso and A bri in (|14|) . Relative to the Wilson operator its eigenvalue spectrum 
is more Ginsparg- Wilson like (cf. Fig. 122]) . and its dispersion relation is more continuum-lik^]. 
As species doubling and global anomalies depend only on topological features of the dispersion 
relation [48,49j, from the conceptual viewpoint this is a Wilson-like fermion. 

When combined with some link smearing and clover improvement, our action was found to 
show good scaling of decay constants even in the physical charm region, and we expect that 
the near-agreement between perturbative and non-perturbative improvement coefficients found 
with the Wilson operator [30],[50l[5T] carries over to this action, too. It appears that lattice 

8 A similar strategy has been adopted for staggered fermions in |47j . 
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Figure 22: Wilson (lapstd,der_std) and Bhllouin (lapJ)ri,der_iso) eigenvalue spectra without 
link smearing and without clover improvement, in 2D (left) and 4D (right). 



perturbation theory is conceptually not any more difficult than for standard Wilson fermions, 
but intermediate expressions may be longer, in particular if several smearing steps are included 
and the backgrounds are made from improved glue [52] . 

Regarding the cost of a simulation with the Brillouin operator it is hard to make generic 
statements. What can be compared is the number of forward applications needed (at a given 
value of M v , cf. Fig.[15l where our Brillouin operator is seen to fare better). However, the cost 
of an individual forward application depends vey much on the architecture used. One extreme 
case is a serial machine which is CPU-limited; in this case 80 neighbor couplings instead of 8 
make each application a factor O(10) slower, whereupon the advantage is gone. On the other 
hand, highly threaded architectures such as GPUs (for an early application to lattice QCD 
see [53]) may be entirely bandwidth-limited; in such a case clever coding might keep the cost of 
a forward application essentially unchanged, relative to the Wilson operator. In our view, the 
upshot is that the usefulness of the Brillouin operator should be tested in phenomenological 
applications where all aspects of a formulation play a role, including the onset of the Symanzik 
scaling regime. In addition, there is a faint possibility that the Brillouin operator might be 
more susceptible to multigrid methods to solve for a given right-hand vector. 

Our Brillouin operator is a specific representative of the class of Dirac operators 

D (x,y) = ^2lfnPn(x-y) + X(x-y) (28) 
where the derivative and the Laplacian are expressed as 

Pn{x-y) = Pl[$B+/J,v-$n-A,tf] 

+ p2 2_j \^x+jl+v,y ~ $x-fi+v,y] 
v 

H~ P3 ^""^ \fix+[i+u+p,y $x— p,+p+p,y] 
v,p 

+ Pi X! [fix+fi+u+p+a,y — S X -p.+P+p+9,y\ (29) 
u,p,cr 
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1.852720547 


240/128 


pi 
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0.136846794 


64/432 


Ai 


-1/2 


-0.060757866 


-8/128 


P2 





0.032077284 


16/432 


A 2 





-0.030036032 


-4/128 


P3 





0.011058131 


4/432 


A 3 





-0.015967620 


-2/128 


Pi 





0.004748991 


1/432 


A 4 





-0.008426812 


-1/128 



Table 2: The pi and A, of the derivative and Laplacian parts of three Wilson-type fermions 
(Wilson's version, the hypercube fermion by Bietenholz et al. JI4J . and our Brillouin operator). 



X(x-y) = A S x>y 

+ Ai ^J[&E+/},y ~ fix-p,,y] 

+ A2 2_}Px+[i-{-i>,y ~ $x-p+i>,y 
p,u 

+ a 3 [<W 



+u+p,y $x— p,+u+p,y] 



p,,v,p 



A4 ^ ] [$x+/l+0+p+(T,y $x— p+u+p+cr,y\ 



(30) 



p,u,p,a 



with the understanding that the sums extend over positive and negative directions mutually 
orthogonal to each other (and in case of the derivative terms also to fi). As discussed in |54j, 
in order to obtain the correct continuum dispersion relation one requires 



2 P i + 12p 2 + 24p 3 + 16p 4 = 1 
A + 8A1 + 24A 2 + 32A 3 + I6A4 = 



(31) 



which all operators in Tab. [2] obey. In addition, our Brillouin (derJso, lap_bri) operator satisfies 

A + 4Ai - I6A3 - I6A4 = 2 

A - 8A 2 + I6A4 = 2 

A - 4Ai + I6A3 - 16A 4 = 2 

A - 8A1 + 24A 2 - 32A 3 + 16A 4 = 2 (32) 

which means that in the weak coupling limit all doublers are lifted by an equal amount, and 

12p 2 + 48p 3 + 48p 4 -l = (33) 

which ensures that the physical branch of the free field dispersion relation E/p has no 0(a 2 ) 
contribution [31], and that the leading cut-off effects in the deviation of the pressure from the 
Stefan-Boltzmann limit are oc 1/A^ 4 [31]. In other words, our Brillouin operator is expected to 
define a discretized version of QCD with decent bulk thermodynamic properties. 

Acknowledgments: We thank Zoltan Fodor for useful correspondence. Computations were 
performed on the high-performance cluster JUROPA at JSC, with resources allocated through 
a NIC grant. Both authors were supported by the SFB TR-55. 
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App. A: Discrete Laplacians and derivatives in 2D 

We give four Laplace stencils in 2D along with their momentum-space representation: 
• Standard Laplacian in 2D: 
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A = 2 cos(fci) + 2 cos(/c 2 ) — 4 [to be read as a 2 A = 2 cos(api) + 2 cos(ap 2 ) — 4] 

• Tilted Laplacian in 2D: 
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A = 2cos(A;i) cos(k 2 ) — 2 [mind the second zero at ki — k 2 — 7r] 

• Brillouin Laplacian in 2D: 
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A = 4cos 2 (A;i/2) cos 2 (A;2/2) — 4 [takes constant value —4 at boundary of BZ] 

• Isotropic Laplacian in 2D: 
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A = [2cos(/ci) cos(k 2 ) + 4cos(A;i) + 4cos(A;2) - 10]/3 [a.k.a. "Mehrstellen" Laplacian] 
We give three x-derivative stencils in 2D along with their momentum-space representation: 
• Standard ^-derivative in 2D: 
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d x = isin(/ci) [to be read as ad x = isin(api)] 

• Brillouin x-derivative in 2D: 
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d x = isin(/c 1 )[cos(/c 2 ) + l]/2 
• Isotropic ^-derivative in 2D: 
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d x = isin(fci)[cos(A;2) + 2]/3 
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App. B: Discrete Laplacians and derivatives in 3D 



We give four Laplace stencils in 3D along with their momentum-space representation: 
• Standard Laplacian in 3D: 
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A = 2 cos(fci) + 2 cos(fc 2 ) + 2 cos^) - 6 
Tilted Laplacian in 3D: 
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A = 2 cos(/ci) cos^) cos^) — 2 
Brillouin Laplacian in 3D: 



[mind the three additional zeros] 




A = 4cos 2 (A; 1 /2) cos 2 (A;2/2) cos 2 (A;3/2) - 4 
Isotropic Laplacian in 3D: 



[takes constant value at boundary of BZ] 
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/48 



A = [cos(fci) cos(k 2 ) cos(k s ) + 3 cos(fci) cos^) + ... + 5 cos(fci) + ... — 25]/6 
We give three x-derivative stencils in 3D along with their momentum-space representation: 
• Standard ^-derivative in 3D: 
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d x = isin(/ci) 

Brillouin ^-derivative in 3D: 




d x = isin(fci)[cos(fc 2 ) + l][cos(fc 3 ) + l]/4 
Isotropic ^-derivative in 3D: 




d x = isin(A;i)[cos(A;2) + 2][cos(A;3) + 2]/9 
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App. C: Discrete Laplacians and derivatives in 4D 

We give four Laplace stencils in 4D along with their momentum-space representation: 



• Standard Laplacian in 4D: 
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• Tilted Laplacian in 4D: 
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A = 2 cos(/ci) cos(A;2) cos^) cos(/c4) — 2 [mind the seven additional zeros] 

• Brillouin Laplacian in 4D: 



1 2 1 

2 4 2 
1 2 1 


2 4 2 
4 8 4 
2 4 2 


1 2 1 

2 4 2 
1 2 1 


2 4 2 
4 8 4 
2 4 2 


4 8 4 
8 -240 8 
4 8 4 


2 4 2 
4 8 4 
2 4 2 


1 2 1 

2 4 2 
1 2 1 


2 4 2 
4 8 4 
2 4 2 


1 2 1 

2 4 2 
1 2 1 



A = 4cos 2 (A;i/2) cos 2 (A;2/2) cos 2 (A; 3 /2) cos 2 (A;4/2) - 4 [constant at boundary of BZ] 



• Isotropic Laplacian in 4D: 



1 


7 


1 


7 


40 


7 


1 


7 


1 


7 


40 


7 


40 


100 


40 


7 


40 


7 


1 


7 


1 


7 


40 


7 


1 


7 


1 
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7 40 7 
40 100 40 
7 40 7 


40 100 40 
100 -2000 100 
40 100 40 


7 40 7 
40 100 40 
7 40 7 


1 7 1 
7 40 7 
1 7 1 


7 40 7 
40 100 40 
7 40 7 


1 7 1 
7 40 7 
1 7 1 



/432 



A = [2cic 2 c 3 c 4 + 7c!C 2 c 3 + ... + 20cic 2 + ... + 25ci + ... - 250]/54 [with ci = cos(/ci) etc. 
We give three x-derivative stencils in 4D along with their momentum-space representation: 
• Standard x-derivative in 4D: 
























-1 1 























4 = isin(A;i) 

Brillouin x-derivative in 4D: 


-1 1 
-2 2 
-1 1 


-2 2 
-4 4 
-2 2 


-1 1 
-2 2 
-1 1 


-2 2 
-4 4 
-2 2 


-4 4 
-8 8 
-4 4 


-2 2 
-4 4 
-2 2 


-1 1 
-2 2 
-1 1 


-2 2 
-4 4 
-2 2 


-1 1 
-2 2 
-1 1 



/2 



/128 



4 = isin(A;i)[cos(A;2) + l][cos(A;3) + l][cos(fc 4 ) + l]/8 
• Isotropic x-derivative in 4D: 



-1 1 
-4 4 
-1 1 


-4 4 
-16 16 
-4 4 


-1 1 
-4 4 
-1 1 


-4 4 
-16 16 
-4 4 


-16 16 
-64 64 
-16 16 


-4 4 
-16 16 
-4 4 


-1 1 
-4 4 
-1 1 


-4 4 
-16 16 
-4 4 


-1 1 
-4 4 
-1 1 



/432 



4 = isin(A; 1 )[cos(A;2) + 2][cos(£;3) + 2][cos(A;4) + 2]/27 
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App. D: Isotropic stencils via Kumar's trick 



A somewhat systematic overview that includes "isotropic" Laplacians in 2D and 3D is presented 
in [21]. In this appendix we review a particularly practical approach for deriving an isotropic 
stencil in any dimension, due to Kumar [20] . 

The standard discretization of the first derivative operator in two dimensions (2D) is 



istd 
'i,3 



1 

2a 



(34) 



From a Taylor expansion [of isin(etfci) =iafei(l— a 2 k\/Q+0{a A )) in momentum space or directly 
in position space] one finds that the standard discrete derivative deviates from the continuum 
derivative through 0(a 2 )-suppressed terms. This may be summarized in the form 



\std 



[1 + ^rdxx){ipx)i,3 



(35) 



where both derivatives on the r.h.s. refer to the continuum. A possible strategy to improve 
rotational symmetry is thus to define the discretized first derivative such that the deviation 
from the continuum behavior xaki is the same in either direction, that is through 

(V>,)S = (l + 7fA)(lk)y ( 36 ) 



where again the operators on the r.h.s. refer to the continuum. The idea by Kumar [20] is to 
factor the bracket and to define the discretized "isotropic" first derivative through 



ISO 



(1 + —dyy){l + —d XX ){i) X ) id 



a 



std 



(37) 



where (135]) has been used in the second step. Moreover, we may replace the second derivative 
in the y direction by its simplest discrete version (i.e. the 1/-2/1 stencil operator), since the 
difference is another a 4 term [of which we did not keep track in (I35H37I) anyway]. This gives 



ISO 

i,3 












-1/2 





-1/2 












-1/2 





1/2 




-1 


1 




1 





-1 




-4 


4 


/12 


-1/2 





1/2 




-1 


1 




1 - 


-1JH 


_i - 4^i-i,j 


-n 


H-i,j-i 


/12 



(3f 



where we have used the stencil notation. Compared to the standard discrete derivative, there is 
a spreading in the transverse direction with a factor |/|/|, respectively. It is easy to generalize 
this procedure to higher dimensions, and the pertinent "isotropic" first derivative operators 
have been given in App. A, B, C for 2D, 3D, 4D, respectively. 

The standard discretization of the second derivative operator in two dimensions (2D) is 



(39) 



and from a Taylor expansion one finds that this is equivalent to 

^2 



\std 
'i,3 



12 xx h3 



(40) 
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where both derivatives on the r.h.s. refer to the continuum. The "isotropic" second derivative 
operator follows by deliberately introducing the same discretization error in the y-direction 



(41) 



and Kumar's trick [20J of factorizing (to the order we are interested in) the continuum expression 



(<••„);•' = (1+tA)(1+-cU0A 



12 



12 



y XX Jl,J 



1 • ^'V) (<•'<.< 



std 



/12 (42) 



yields (the simplest) "isotropic" second derivative operator. Compared to the standard discrete 
second derivative in 2D, there is a spreading in the transverse direction by a factor j^/jf/^- 
It is easy to generalize this procedure to higher dimensions, and we shall just give the result: 

• Isotropic second x-derivative in 2D: 



1 


-2 


1 




10 


-20 


10 


/12 


1 


-2 


1 





d 2 



cos(fci) - l][cos(Jfe 2 ) + 5]/3 
Isotropic second x-derivative in 3D: 



1 


-2 1 


10 -20 10 


1 -2 


1 




10 


-20 10 


100-200100 


10 -20 


10 


/144 


1 


-2 1 


10 -20 10 


1 -2 


1 




d 2 = 

^x 


[cos(fci) 


- l][cos(fc 2 ) + 


5][cos(fc 3 ) 


+ 5]/18 



Isotropic second x-derivative in 4D: 



1 -2 1 
10 -20 10 
1 -2 1 


10 -20 10 
100 -200 100 
10 -20 10 


1 -2 1 
10 -20 10 
1 -2 1 


10 -20 10 
100 -200 100 
10 -20 10 


100 -200 100 
1000-20001000 
100 -200 100 


10 -20 10 
100 -200 100 
10 -20 10 


1 -2 1 
10 -20 10 
1 -2 1 


10 -20 10 
100 -200 100 
10 -20 10 


1 -2 1 
10 -20 10 
1 -2 1 



/1728 



d 2 x = [cos(Jfei) - l][cos(fc 2 ) + 5][cos(A;3) + 5][cos(fc 4 ) + 5]/108 

Upon adding the "isotropic" second derivative in the y (and possibly z, t) direction, one gets the 
"isotropic" Laplacian stencil in 2D, 3D, 4D, as given in previous appendices. We emphasize that 
the "isotropic" Laplacian establishes better rotational symmetry near the center of the Brillouin 
zone. What proves most useful in many applications (including our goal of designing more 
continuum-like lattice Dirac operators), however, is isotropy at the boundary of the Brillouin 
zone, and this is achieved through the "Brillouin" Laplacian, as given in App.A,B,C. 
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